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ABSTRACT 

We present a hybrid code combining the OpenMP-parallel tree code VINE with an 
algorithmic chain regularization scheme. The new code, called “rVINE”, aims to signif¬ 
icantly improve the accuracy of close encounters of massive bodies with supermassive 
black holes in galaxy-scale numerical simulations. We demonstrate the capabilities of 
the code by studying two test problems, the sinking of a single massive black hole 
to the centre of a gas-free galaxy due to dynamical friction and the hardening of a 
supermassive black hole binary due to close stellar encounters. We show that results 
obtained with rVINE compare well with NBODY7 for problems with particle numbers 
that can be simulated with NBODY7. In particular, in both NBODY7 and rVINE we 
find a clear N-dependence of the binary hardening rate, a low binary eccentricity and 
moderate eccentricity evolution, as well as the conversion of the galaxy’s inner density 
profile from a cusp to a a core via the ejection of stars at high velocity. The much 
larger number of particles that can be handled by rVINE will open up exciting oppor¬ 
tunities to model stellar dynamics close to SMBHs much more accurately in a realistic 
galactic context. This will help to remedy the inherent limitations of commonly used 
tree solvers to follow the correct dynamical evolution of black holes in galaxy scale 
simulations. 

Key words: black hole physics — stars: kinematics and dynamics — galaxies: evo¬ 
lution — galaxies: nuclei — methods: numerical 


1 INTRODUCTION 


The p resence of supermassive black holes (SMBHs; e.g. iReesI 
1198411 with masses of 10® Mq < Mbh < 10^® Mq hosted in 
the central regions of virtually all massive spheroids in the 
nearby Univ erse, including the Galactic bulge, is now firml y 
established (iRichstone et al.lll998l : iKormenA fc Hdl2013ll . 
SMBHs must have also been present at much earlier phases 
of our Universe, powering active galactic nuclei (AGNs) and 
quasars from only a few hundred milli on years after the 


Big Bang throughout cosmic history (e.g. LyndemBelj ljjGgl : 


iFan et 7 uI 2001| : Tcivano et al.l[2^01ll : iMortlock et al.ll20li l. 
We now think of SMBHs as integral components of 


galactic nuclei, possibly playing a decisive role in shap¬ 
ing the structure and morphology, as well as the gas 
and thus stellar content of massive galaxies. The ACDM 
paradigm for structure formation together with a number 
of surprisingly tight relations between the SMBH masses 
and fundamental properties of the galactic bulges hosting 
them - e.g. the s pheroid luminosity llKor mendv fc Richs tond 
I 1995 II and mass jMagorrian et al.lll99^:lHaring fc Rixl2004ll . 
and the stellar velocity dispersion llGebhardt et al.l 2000l : 
iFerrarese fc MerrittI l2000l : I Tremaine et al.l 2002ll - mggests 
co-evolution of the hier archically growing galaxies and their 
central black holes (e.g. L ynden- Bell 1969 : ISilk fc Reeslll998l : 


iKauffmann fc Haehnelt 20001 ). 
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At present, variants of p article-based Smoothed Parti- 
cle Hydrodynamics (SPH, e.g. lWadslev et al'1l2004l : ISpringell 
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200^, Krid-bas ed A daptive M esh Refinement (AMR, e.g. 


Kravtsov et al.lll997l: lTevssieill200^ 1 codes, or moving mesh 


codes l SDringelll201Cll ). are the methods of choice for numer¬ 
ical simulations of cosmological galaxy formation, exploit¬ 
ing the high dynamic range and spatial flexibility in reso¬ 
lution. The gravity solvers in these codes either employ a 
particle-mesh scheme or are tree based and assume the sim¬ 
ulated system to be collisionless, with two-body relaxation 
time-sc ales exceeding the ag e of the s ystem. In ’tree’ algo- 
rithms dBarnes HutI 119861 . see also iMcMillan fc AarsethI 
I 1993 I for a collisional tree code) the gravitational force on a 
single particle from a distant group of particles is approx¬ 
imated by a multipole expansion about the group’s centre 
of mass, and the N-body particles represent massive tracer 
particles that sample the underlying smooth gravitational 
potentials. To reduce the graininess of the potential, gravi¬ 
tational forces are ’softened’ on small spatial scales and the 
softening length e is the natural resolution limit of the code. 
Unfortunately, this also means that close two-body encoun¬ 
ters with a massive body such as a supermassive black hole 
can - by construction - not be computed accurately. An alter¬ 
native, much more cost-intensive approach of calculating the 
gravitational forces in an astrophysical system is the direct 
summation of each particle’s gravity on every other particle 
in the system. Combined with high-order integrators, this 
method is a very accurate way to calculate the gravitational 
forces and is widely used to simulate collisional N -body sys¬ 
tems (e.g. lAarsethlllg^. l2003bl: iHut et al. IH). 

Owing to the inherent limitations in the numerical 
methods, studies of stellar dynamics in the vicinity of 
SMBHs to date could only probe single separate aspects 
of the full problem. On the one hand, direct N-body simu- 


galaxies or merger remnants (e.g. lEbisuzaki et al.l 

199l|; [7l 

Milosavlievic & Merritt! 200ll. 20031: 

Berczik et al.l 

20o4 

Merritt et al.l 20071: Khan et al. 2011 

. 2 OI 2 I: Preto et al.l 

2 OIIII use idealized initial conditions to represent the in- 


ner parts of galaxies in which the massive binaries are then 
embedded and evolved. Due to the steep scaling of required 
computing time with particle number of order 0{N‘^) these 
studies are still rather limited in the particle number, hin¬ 
dering a self-consistent treatment of full galactic environ¬ 
ments. On the other hand, simulations of galaxy mergers 
or cosmological simulations of structure formation includ¬ 
ing SMBHs (e.g.J Springel et al.||2 005l: iDi ^Matteo_e^^ ]|2008l: 


Johansson et al. 120091: Booth fc Schaydl2 009l : Martizzi et al.l 

2 OI 2 I : IChoi et al.ll2012l : ISiiacki et al.ll20l4 l have difficulties 


to capture the dynamics of SMBHs and their surrounding 
stars below the resolution limit. This leads to uncertainties, 
e.g. in the dynamical friction time-scales for the SMBHs, 
and affects the density- and velocity-profiles of the stellar 
background through interactions with the SMBHs, and the 
hardening and merging time-scales of close SMBH pairs. 
The latter is generally assumed ’a priori’ to happen fast in 
these simulations, much reducing the accuracy and predic¬ 
tive power of such simulations. Hence, there is still substan¬ 
tial uncertainty in the current understanding of the dynam¬ 
ical evolution of SMBHs and their surrounding star clusters 
in realistic cosmological settings, which directly feeds back 
into uncertainties in our understanding of how SMBH sin¬ 
gles or multiples influence the structure of galaxies. 

The main goal of this paper is to help remedy these 


shortcomings by combining the best parts of the two numeri¬ 
cal approaches: a regularization method to efficiently and ac¬ 
curately compute the dynamics close to the black holes and 
a fast tree code to treat the global galactic dynamics. This 
goes in line with the development of similar recent hybrid 
codes (as discussed in the next Section) and a software in¬ 
terface designed to efficiently combine different stan d-alone 
code architectures dPortegies Zwart et al.|[^009l . l2013tl . With 
the new algorithm, we will be in a position to better take 
into account the relevant dynamical processes regarding the 
interaction between SMBHs and the stars in their environ¬ 
ments and other (SM)BHs, in principal without limitations 
on the spatial and temporal resolution down to scales where 
other types of physical phenomena become important, e.g. 
gravitational wave induced coalescence of SMBH binaries 
or the tidal disruption of low-angul ar-momentum stars (e.g. 
|Pretoriu3 l2005l : iLodato et al.l l2009l l . This might be an im¬ 
portant next step towards investigating the dynamical co¬ 
evolution of SMBHs and their host galaxy nuclei in a self- 
consistent manner in galaxy-scale or cosmological simula¬ 
tions. 

In this paper, we present the details of our hybrid 
N-body code which will help us to focus on the role of 
gravitational dynamics in the interplay between supermas¬ 
sive black holes and realistic representations of the central 
regions of their host galaxies. The paper is structured as 
follows. In Section [2] we describe the details and structure 
of the new hybrid code “rVINE” and test its performance 
in Section |4] after we have described the numerical set-up 
in Section [31 First tests on the code in comparison with the 
direct N-body code NBODY7 and the tree code VINE are 
presented in Section [5] We discuss our results in Section [6] 
and, finally, summarize and draw our conclusions in Section 


2 A NEW REGULARIZED TREE CODE 

In this Section we present the structure of the regular¬ 
ized tree code called “rVINE”. A number of currently avail¬ 
able N-body codes employ regularization techniques in¬ 
tended for the integration of strong gravitational interac¬ 
tions but are primarily developed for integrations of col- 


tral reEions of Ealactic nuclei (see AarsethI 2003a. 2007: 

Harfst et al.l 20081: Gaburov et al.l 20091: 

Nitadori & AarsethI 

2 OI 2 I: Berczik et al. 20131: IWans et al.l 

2015f). In addition. 


there are two recent hybrid codes combining t ree a nd N- 
body codes. The BRIDGE code (iFuiii et al.l l2007fl uses 
a sim ple fixed time-step oct-tree, while the Bonsai tree 
code (iBedorf et al.l 1201^ is an oct-tree run entirely on 
GPUs. Both hybrids use a symplectic mapping method 
ll Wisdom fc Holmanlll99ih to couple the tree to a direct N- 
body algorithm with a fourth-order Hermite integrator. 

In rVINE we compute the evolution of a sub-system 
of particles near the black hole by means of a regulariza¬ 
tion method, while regions of the galaxy further out with 
long relaxation times and basically unaffected by the pres¬ 
ence of the black hole, are integrated by a collisionless tree 
code. To this end, we combine two publi shed algorithms: 
a version of the tree/SPH code VINE llWetzstein et al.l 
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Figure 1. Projected trajectories of regularized particles in the chain around a supermassive black hole. The trajectories are shown in a 
time interval At 1, which corresponds to one crossing time of the subsystem shown. 


I 2 OQ 9 I : iNelson et al.l l20Q9h and the algorithmic re gulariza¬ 
tion (AR) chain method ([Mikkola &: Aarsethll2002l . see also 
iHellstrom &: Mikkolall20ig for a general discussion). The lat¬ 
ter was kindly provided as a stand-alone code by Seppo 
Mikkola. 

VINE is an OpenMP-parallelized tree/SPH code em¬ 
ploying a binary tree algorithm and an individual hi¬ 
erarchical block t ime-step scheme (jWetzstein et al.l l2Q09l : 
iNelson et al.l l2009lf1 . The AR-chain method is an efficient 


and extremely accurate method to study close dynamical 
few-body encounters and is capable of handling even (re¬ 
peat edj_two 2 body_conisions jMikkola^^_^ajse^ 200i see 


also IPreto &: Tremain3l 19991 : iMikkola fc Tanikawal 1999al f9l . 

This is achieved by effectively removing any singular be¬ 
haviour in the equations of motion by a time transfor- 
mation in the Hamiltoni an of the regularized sub-system 
ijMikkola fc Aarsethll20o3 FL The coordinates and the (orig¬ 
inal) time and their respective ’momenta’ are integrated us¬ 
ing a simple leapfrog method. In addition, a Bulirsch-Stoer 
extrapolation method l|Gragelll96lTl : iBulirsch fc Stoeilll96fil l 
is applied to guarantee high accuracy, as well as a chain 
concep t of smallest inter-particle vecto rs to reduce round-off 
errors llMikkola fc Aarsethlll990l . Il993l f . In our present ver¬ 
sion of the AR-chain, chain particles are sorted according to 
their gravitational forces, not according to their distance. It 
also includes a method to handle velocity-dependent forces, 
which allows us, in principle, to treat additional viscous 
and relativistic terms in the regularized force calculations 
llMikkola fc Merritdbood l. The new code, however, is purely 
Newtonian at the present stage. 

There exist other regularization schemes which are com¬ 
parable to the AR-chain in accuracy a nd speed, e.g. the 
KS-whee lspoke j[Zax^jjj7^jAaigetlJ l200^ 1 and the KS-chain 
method llMikkola fc Aarsethl Il99^ . Due to limitations of 


^ Note that in the present paper we will discuss new develop¬ 
ments done in the parts of VINE related with the leapfrog inte¬ 
grator, an individual time-step scheme and no Smoothed Particle 
Hydrodynamics. 

^ However, the singularities do formally remain in the trans- 
formed equations of motion - unlike in schemes based on 
iKustaanheimo fc Stiefell ll 19651 . (KS)) regularization, which ap¬ 
plies a time and a coordinate transformation. 


these methods in the context of stellar dynamics around 
one or several SMBHs we decided to use the AR-chain as 
our principal regularization algorithm. For example, large 
mass ratios may lead to a loss in numerical accuracy for the 
less massive bodies in a KS-chain, whereas the wheelspoke 
has difficulties in treating multiple heavy bodies on an equal 
footing. 

The gravitational forces for the majority of the particles 
are computed with VINE’S fast binary tree scheme with¬ 
out regularization, using a pre-defined spline or Plummer 
softening, while particles near a designated massive particle 
(SMBH) become members of a compact subsystem which 
is integrated in the AR-chain in its centre-of-mass reference 
frame. 

Chain integration starts if any particle comes closer to the 
SMBH than rj^BH < rdiain.o, where rchain.o defines the initial 
chain size and is an input parameter which has to be chosen 
at the start of a simulation as described in the following. 

The intended purpose for including particles in the 
chain integration is twofold. Firstly, we want to accurately 
follow close orbits near the black holes, i.e. within a fair 
fraction of the SMBHs’ influence radii, and secondly, we 
need to overcome the limitations posed by the gravita¬ 
tional softening, i.e. for encounters within ~ a few times 
e = max(£BH,e*), where e* and £bh denote the gravita¬ 
tional softening lengths of the stellar particles and SMBHs, 
respectively. Hence, we determine the initial chain radius, 

fchain.o = max(Q ■ Tinfl, P ■ s), (1) 

at the start of the chain, where a and /3 are input parameters 
and Tinfi is the gravitational influence radius of the SMBH. 
For practical purposes, we define rinfl = min(rf,J^, using 
two commonly used proxies for the gravitational influence 
radius of the SMBH, being 1) the radius enclosing twice 
the mass of the SMBH, = r(< 2 Mbh), and 2) the 
radius within which the gravitational force of the SMBH 
dominates over the self-gravity of the stellar background, 
D'nfl — GMbh/u^. Here, Mbh is the mass of the black hole 
and G the gravitational constant. The velocity dispersion 
cr is determined by averaging over the nearest 50 particles 
outside the chain. We generally find good results for 1 ^ 
a ^ 1.5 and 1 ^ /9 ^ 2. 

The chain’s center-of-mass is treated as a massive parti- 
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Figure 2. Illustration of the different integration regions near a regularized massive particle in the hybrid AR-regularized tree code. The 
centre-of-mass reference body, consisting of the members within the regularized subsystem (r < rchain)i is surrounded by a swarm of 
nearby particles (’perturbers ’) which are considered as external force terms in the chain force calculations (r < rpert, where the critical 
radius rpert depends on a tidal criterion) and, themselves, experience direct N-body forces from the resolved chain (see text). Further 
out, we indicate the regime where the direct integration of particles (open circles) switches to the multipole approximations, depending 
on the acceptance criteria in the tree code (indicated by the groups of particles in the large star symbols). 


cle in the tre43, i.e. it is included in the tree force calculations 
and advanced in time within the tree code. The chain parti¬ 
cle is advanced in time on the smallest tree time-step and we 
formally set the gravitational softening of the chain particle 
to zero in the tree-code. Tree particles that become mem¬ 
bers of the chain are converted into ’ghost’ particles with no 
further advancement in the tree. This is done by assigning a 
very small, but finite-sized mass in the tree code data struc¬ 
ture, rendering their contribution to the gravitational forces 
negligible. Furthermore, ghost particles are not allowed to 
determine the size of the time-step in the tree. If the chain is 
active, the member particles within the chain are advanced 
using the AR-chain integration, every time particles in the 
tree code on the smallest time-step level are being advanced. 

The equations of motion of the chain members include 
external forces exerted by a set of nearby tree particles we 
call “perturbers”, which are identified via a tidal criterion, 


rj.CoM < 


/ 2 rrij \ 3 

( "FT 1 ^ ^’ 

VTcrit fWchain / 


( 2 ) 


where rrij is the mass of perturber j, Mchain is the total mass 
in the chain, and 7crit a dimensionless parameter which de¬ 
fines the relative tidal perturbation on the chain. We typi¬ 
cally set rcrit = min(rchain,rchain,o) in the simulations pre¬ 
sented here. For better accuracy, the perturber positions 
relative to the center-of-mass particle are predicted (to 1st 
order) to the current time at each force calculation in the 
AR-chain. The perturber forces have to be predicted many 
times during the numerous sub-step cycles of the Bulirsch- 
Stoer extrapolations. To improve the overall performance of 
the chain part of the code we, therefore, have implemented 
parallel routines for the predictions of the perturbers and 
the computation of the perturber forces. 


® Henceforth we will call this particle simply the “chain particle”, 
denoting the chain’s centre-of-mass particle that is advanced in 
the tree code; not to be confused with a single “particle in the 
chain”, which we will equally call a “chain member” from now on. 


Table 2. Parameters for the rVINE calculations. 


Simulation^ a 

/3 

7 

'Tcrit 

A rVine 1.0 

1.0 

1.5 

10-^ 

B rVine 0.1 — 0.3 

1.0 - 4.0 

1.5 

10-4/10-® 


“ Note that, throughout the text, added suffixes will state the actual 
particle numbers of the simulations. 


Likewise, the force contributions of individual (’re¬ 
solved’) chain particles are added to the gravitational forces 
for the perturber particles during the force updates in the 
tree. The gravitational force on the chain particle is also cor¬ 
rected by resolving the chain particles. If the chain particle 
is active, we first calculate the gravitational force exerted 
on it by particles within the tree code, but subtract the tree 
forces (i.e. direct, softened N-body interactions) from the 
perturber particles again. Then the (direct) gravitational 
force from the perturbers on each individual chain particle 
Is calculated, and added up as a correction to the chain par¬ 
ticle’s force according to 


.^chain 


C^CoMjCorr — 


Md 


- E 


^^chain,i Cf-chain,i* 


(3) 


After advancing the chain one full time-step its member¬ 
ship is updated. Perturber particles are added to the chain 
if they come closer than the “chain radius”, i.e. 


^j,CoM ^ Mchain; 


(4) 


where we define the chain radius as the largest distance of 
any chain member relative to the chain centre-of-mass in the 
last chain step. Particles are removed from the AR-chain if 
they recede far enough from the chain’s center-of-mass 


^j,CoM > 7 X rdiain.O — ^escape 


(5) 
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Table 1. Parameters for the different N-body calculations. 


SimulatioiA’^ 

w 

Mtot 

VQ 

mi,/M 


Abh 

^init 

'bh 

..init 

AfBH/Aftot 


A NbodyT 

look 

1.0 

1.0 

10-5 

- 

1 

2.14 

0.46 

10-5 

- 

A Vine-El/A Gadget-El 

100k 

1.0 

1.0 

10-5 

0.02 

1 

2.14 

0.46 

10-5 

0.1 

A Vine-E2/A Gadget-E2 

100k 

1.0 

1.0 

10-5 

0.02 

1 

2.14 

0.46 

10-5 

0.02 

A rVine 

100k 

1.0 

1.0 

10-5 

0.02 

1 

2.14 

0.46 

10-5 

- 

B Nbody7 

lOk-lOOk 

1.0 

1.0 

10-"^ - 10-5 

- 

2 

0.1/0 

0.28/0 

5 X 10-5 

- 

B Vine 

lOk-lM 

1.0 

1.0 

1 

o 

1 

1 

o 

0.01 

2 

0.1/0 

0.28/0 

5 X 10-5 

0.01 

B rVine 

lOk-lM 

1.0 

1.0 

1 

O 

1 

1 

o 

0.01 

2 

0.1/0 

0.28/0 

5 X 10-5 

0.01 


“ All quantities are given in code units. 

^ Note that, throughout the text, added suffixes will state the actual particle numbers of the simulations. 



Figure 3. Projections of the initial stellar surface densities (in code units) for a realisation of model A_rVine_100k. The orbital evolution 
of a SMBH, placed on a circular orbit at the half-mass radius, is shown as the black solid line and symbols for the first orbital period. 


and move radially away from the chain centre. We typi¬ 
cally choose 7 = 1.5. Upon removal of a chain member, its 
mass, position and velocities in the global (i.e. tree) reference 
frame are restored. If the number of chain members would 
fall below A);hain < 2 after an escape the chain is terminated 
and the remaining two chain members are restored in the 
tree. Likewise, the chain is terminated in the (unlikely) case 
that the designated massive particle is removed from the 
chain — unless the chain membership includes several mas¬ 
sive particles, in which case one of the remaining SMBHs 
is chosen to be the new centre-of-mass particle in the tree- 
code. A representative example of near-Keplerian, perturbed 
orbits of particles in the chain is illustrated in Figure [T] 

The regularization of particle orbits near a designated 
(massive) particle, leads to different integration regions for 
different particles in the code, as illustrated in Figure [2] 
With increasing distance to the chain’s centre-of-mass the 
particles are either 

(i) regularized members of the chain (within the dashed 
circle in Figure [2|, feeling external gravitational forces from 
the perturbers (filled circles) only, 

(ii) perturber particles which feel the gravitational forces 
of the individual chain members, or 

(iii) tree particles, whose gravitational forces onto the 
chain particle may either be calculated via direct summation 
(open circles) or approximated by a multipole expansion of 


a group of particles (indicated by star symbols), depending 
on the values for the acceptance criteria chosen in the tree. 

For cost reasons, we restrict the total number of chain 
members and perturbers. We have tested up to values of 
~ 250 and = 5000 without any loss in the sta¬ 

bility of the code. 

In short, a typical time-step in the hybrid code proceeds 
as follows. 

(1) Determine the next time-step and active particles for 
the tree. 

(2) If the chain is not active, or, the system does not evolve 
on the smallest time-step, continue with step (6). 

(3) Integrate the members of the chain using the AR-chain 
method. Treat the external forces exerted by nearby 
particles (“perturbers”) as perturbations to the members’ 
gravitational forces. 

(4) Check for absorption by and escape from the chain. 
Perform a search to identify the perturbers via a tidal 
criterion in regular time intervals (see Equation [^. 

(5) Terminate the chain if Achain < 2. 

(6) Perform regular leapfrog time-step in the tree code in 
the ’drift-kick-drift’ (DKD) scheme: 

(a) Update tree particle positions at the half time-step. 

(b) Compute gravitational forces for the active tree parti¬ 
cles. If the chain particle or any perturber is active include 
force corrections due to the gravitational forces between the 


© 2012 RAS, MNRAS OOO.IllfTTl 


















6 S. J. Karl, S. J. Aarseth, T. Naab, M. G. Haehnelt, R. Spurzem 


perturber particles and the resolved members of the chain, 
(c) Update tree particle positions and velocities at the full 
time-step. 

(7) Update tree particle time-steps in the individual time- 
step scheme and update the tree structure if necessary. 

(8) If the chain is not active check for particles near the 
SMBH that fulfill the conditions to begin chain regulariza¬ 
tion. 


3 NUMERICAL SET-UP 

In this Section we shortly describe the numerical set-up and 
the different numerical models we will introduce in the fol¬ 
lowing Sections. In particular, we will detail the way the 
code works by discussing an example simulation and testing 
its performance for a number of different code parameters 
in SectionlD In Section[5]we test the code against a number 
of other currently available N-body codes, such as the VINE 
and Gadget-3 tree codes and the NBODY7 direct summa¬ 
tion code. All of the rVINE, VINE, or Gadget-3 simulations 
presented in this paper were run on the COSMOS cluster at 
DAMTP, Cambridge. For the NBODY7 simulations we used 
single nodes on the Wilkes cluster at the High Performance 
Computing (HPC) Service of the University of Cambridge, 
consisting of a Dell PowerEdge T620 server a 12 Intel Xeon 
E5-2670 CPUs plus two NVIDIA K20 GPUs. 

As ou r princ ipal numerical model we use a non-rotating 
iHernauis^ lll99d l sphere to represent the galactic nucleus. 
The Hernquist density profile follows a p{r) oc r~^ power- 
law at small radii (r << ro) while converging quickly, with 
p(r) oc r~‘^, to a finite mass at large radii (r >> ro). All 
models are set-up with unit total mass (Mtot = 1) and 
scale radius (ro = 1), and the gravitational constant G is 
set to unity in all codes used throughout this paper. Note 
that, for convenience, we will use a system of code units 
in the plots shown throughout the paper. Identifying our 
model, for instance, with a small spherical (dE/SO) galaxy 
or a galactic nucleus of mass Mtot = 10^° Mq and scale 
length ro = 1 kpc yields time units and velocity units of 
~ 4.7 Myr and ~ 207kms“^, respectively. In code units 
the half-mass radius of the Hernquist sphere is given as 
’’' 1/2 = (l + t/S) ro ~ 2.41 with a half-mass dynamical time of 
fdyn.h = 27. In the different simulations considered here the 
galaxy is realised with total particle numbers in the range 
10 "* — 10® of equal-mass stellar particles (m* oc N~^) for 
VINE and rVINE0 For NBODY7, however, we restrict our¬ 
selves to N ^ 10® particles due to the steep 0{N^) scal¬ 
ing of the direct N-body code. Additionally, we place a few 
massive particles representing the SMBHs in the simulation, 
with initial conditions depending on the specific simulation 
run. Models using a single SMBH are denoted as models 
’A’ (Sections U] and 15.111 while models using two SMBHs 
are denoted as models ’B’ fSections l4l and 15.211 . To reduce 
the scatter in our simulations, we will show the results for 
all models as averages over several independent Monte-Carlo 
realisations of initial conditions with different random seeds. 

^ Note also that we plan to employ much higher particle numbers 
in rVlNE, well above what we have used here for the comparison 
tests, in future (see Section [bjl. 



Figure 4. Time evolution of the number of members in the 
chain and the number of perturbers in a realisation of model 
A_rVine_100k, with one SMBH initially set on a circular orbit 
(see Section IS.lll . Average numbers are given as text and indicated 
by the dashed lines. 


Depending on V, we set up four different realisations for 20k 
^ V ^ 100k, two realisations for 100k < N < IM and one 
realisation for N — IM. Details of the set-up together with 
the explicit name for each simulation run are given in Table 
[T] As an example, we show the early orbital evolution of a 
SMBH on the underlying initial stellar surface density for a 
particular realisation of models ’A_rVine_100k’ in Figure 
[3] (see also Section 15.11) . 

In all simulations using either VINE or rVINE , we intro¬ 
duce a simple Plummer softening llAarse thlll963l ) with dif¬ 
ferent values for the star particles e* and the SMBHs (ebh). 
However, interactions between chain members as well as in¬ 
teractions with the (active) chain centre-of-mass particle in 
the tree code are not softened at all. For the force calcu¬ 
lations in the tree we use a “Gadget” multipole acceptance 
criterion (MAC) with error tolerance 6 — 10~® and conserva¬ 
tive values for the accuracy parameters race,* = rvei,* = 0.3 
and Taec.BH = TVei.BH = 0.003 for the leapfrog integration 
fsee lWetzstein et al.ll200^ . In the Gadget-3 simulations the 
same integration accuracy, r = 0.02, is adopted for both 
stellar and BH particles, along with an error tolerance of 
5x 10~®. For NBODY7 we set the chain radius to 1.25 x 10“® 
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Figure 5. Time evolution of characteristic radii in a realisation of 
model A_rVine_100k, with one SMBH set on a circular orbit (see 
Section ISTT I: the radius of the chain (red solid line), the distance 
of the nearest perturber to the chain centre-of-mass (green solid 
line) and the radius of escaping chain particles (black dashed line). 

and do not make use of softened gravitational interactions. 
The total accumulated relative energy error stayed below 
< 10~^ for all simulations with an accuracy parameter of 
To = 0.02. The parameter ranges used to define the per¬ 
turber particles and chain members (see Section mi in the 
rVINE simulations can be found in Table [2]. 


4 TESTING THE CODE AND CODE 

PERFORMANCE 

In this section, we carry out some basic examples and tests 
on the code performance of the AR-regularized tree code. 

In Figure |4] we illustrate the time evolution of the 
number of particles in the chain (upper panel) and the 
number of perturbers (lower panel) for one realisation of 
model A_rVine_100k, following the orbital evolution of a 
SMBH in a Hernquist sphere consisting of A = 10® parti¬ 
cles (see also Section ICT) . The SMBH is assigned a mass 
of Mbh = 10“® Mtot and set on an initially circular orbit 
at the half-mass radius (rBH(t = 0) = 2.41). As the SMBH 
sinks to the centre the number of chain members as well 
as perturbers increases significantly at time t ~ 520, i.e. at 
about the time the SMBH sinks to the central high-density 
regions within < 10% of the half-mass radius. Shortly after 
we reach the maximum number of 55 particles in the chain 
and a maximum number of ~ 780 perturbers, while the av¬ 
erage numbers are only (Achain) ~ 20 and (Apert) ~ 500, 
respectively. This highlights the need to perform some ex¬ 
ploratory simulations to determine whether rVINE can actu¬ 
ally handle the range of particle densities around the SMBH 
encountered over the simulated time span. 

Overall, a fraction of 6.9% of the 10® particles, was sub¬ 
ject to integration in the chain during this run. Chain inte¬ 
gration was active for a total of 42% of the total run time, 
tmax = 800. At the early stages of the simulation, the chain 


is used only intermittently to treat the occasional strong 
stellar encounters near the SMBH in the low density envi¬ 
ronment. At later phases, when the SMBH is near to the 
centre, the chain is used more intensively, with the chain 
being active continuously for ~ 285 time units. Stellar par¬ 
ticles are typically included repeatedly in the chain with an 
average number of ~ 8 recurrences. 

Figure [S] shows the corresponding time evolution of 
some characteristic radii from the simulation shown in Fig¬ 
ure [4] The escape radius (black dashed line), given by Equa¬ 
tion serves as an effective upper bound for the chain 
radius (red solid line). For times t < 300 the minimum dis¬ 
turber distance, rp,miii, (green solid line) is generally well 
above the escape radius, before some perturbers may come 
closer to a chain that has only a few and, by chance, quite 
compact members while the SMBH is sinking to denser cen¬ 
tral regions (300 < t < 520). Once the SMBH is close to the 
centre of the Hernquist sphere (t > 550) and the number 
of particles in the chain has increased by a factor of ten, 
f'chain and Tp^min both oscillate around Teacape, which nat¬ 
urally arises when a number of particles lives close to the 
conditions for both absorption and escape being satisfied 
at a certain time. In this case, in rVINE we prioritise the 
absorption of near perturbers over a (delayed) escape of a 
chain member, until the chain radius has grown by 5% over 
the nominal escape radius. A noticeable oscillation around 
the escape radius can then occur if a series of subsequent 
absorptions of perturbers near the SMBH takes place. 

In Figures |S] and [7] we investigate the performance char¬ 
acteristics of rVINE using a set of simulations of our basic 
Hernquist model with two SMBHs, each having a mass of 
Mbh = 5 X 10“® Mtot. One of the SMBHs is initially set 
on a circular orbit close to the centre (tbhi = ±0.1 and 
Uy,BHi ~ 0.28), the other one is at rest at the origin (mod¬ 
els B_Vine and B rVine; see also Section EH]). All runs 
shown in Figures |6] and [7] were performed on a single node 
(8 CPUs) of Cosmos2, except for a few comparison runs 
using NBODY7 (see Figure El upper panel) which were per¬ 
formed using 8 CPUs plus acceleration from two NVIDIA 
K20 GPUs on the Wilkes cluster. 

The three panels in Figure El show the wall clock time 
required to evolve the simulation to t = 1 (upper panel), as 
well as the average number of particles in the chain (middle 
panel) and the average number of perturber particles (bot¬ 
tom panel) as a function of the initial chain radius Tchain.o 
and particle number A. Interestingly, the run time does not 
seem to depend strongly on the choice of the initial chain 
radius; it typically changes by a factor of ~ a few, and at 
most by a factor of ~ 8, when varying Tchain.o up to a fac¬ 
tor of 3. In addition, the scaling with A is still relatively 
shallow for all initial chain radii as the computing time is 
dominated by the tree. For comparison, we show the scaling 
obtained with NBODY7 using the same initial conditions 
for A ^ 2 X 10® (black dashed line). Due to the fact that 
we can use the additional acceleration of the two NVIDIA 
K20 GPUs (plus some contribution from running on a dif¬ 
ferent system) the total computing time is actually a factor 
of ~ 5 to ~ 10 times lower compared to rVINE at the lowest 
particle numbers. However, the steeper scaling of NBODY7 
(Twaii oc N^'^) yields comparable computing times already 
for A > 2 X 10®. Extrapolating this scaling would give clear 
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Figure 7. Performance characteristics for the rVINE code as a 
function of varying initial chain radii for the model used in Figure 
Shown are the code run time (top panel) and the relative energy 
error (bottom panel) versus the accuracy parameter of the tree 
code time-step criteria, r*. All runs were evolved for one code 
time unit (At = 1). 



Figure 6. Performance characteristics for the rVINE code as 
a function of varying initial chain radii in calculations of two 
SMBHs at the centre of a Hernquist sphere, one being set on a 
circular orbit at o^bhi = 0.1 and the other one being at rest at 
the origin (see also Section lot . The three panels show the code 
run time in seconds {top), the average number of particles in the 
chain (middle), and the average number of perturbers (bottom) 
versus the total number of stellar particles in the simulations. All 
runs were evolved for one code time unit (At = 1). 


advantages to rVINE for N > 10® particles in that particular 
case. 

The average number of particles in the chain shows a 
roughly linear scaling with total particle number for N > 2x 
10 “* while the average number of perturbers has a shallower 
N-dependence. The latter can be understood by considering 
the fact that the region of the perturber particles actually 
becomes smaller with increasing N, since the mass of the 
perturbers scales as oc N~^ (see Equation (HJ). Hence, if the 
number of particles in the chain were to scale strictly linearly 
with N, or in the limit of the SMBH dominating the total 
mass of the chain, e.g. for a high SMBH-to-star mass ratio 
and small rchain.o, we would expect (Vpert) to be largely 
independent of N. The observed shallow increase of (Vpert) 
with N is likely to be caused by the non-trivial non-llnear 
scaling (Vchain) observed in the middle panel. Both (Vchain) 
and (Vpert) show a scaling going roughly as V oc r^i,ain, 0 j 
as expected for the c entral parts of t he Hernquist profile, 
where M(r) oc (see lHernauistlll990l . their equation (3)). 

The upper panel of Figure [7] shows the total run time of 
the simulations as a function of the accuracy parameter for 
the tree code time-step criteria, r*. In principle one is free 
to choose different values for the different time-step criteria 
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used in VINE fsee S ection [3] and equations (10) - (12) in 
IWetzstein et al.ll2009li . but we decided to adopt one, identi¬ 
cal value of T* for all criteria and fixed the particle number at 
N — 10®. The time-step size increases linearly with r* lead¬ 
ing to an overall decrease in the run time for all values of 
rchain.o, albeit with some non-negligible scatter. Within the 
scatter, the total wall time scales roughly as T^aii oc Tchain.o 
for fixed t*. In the lower panel of Figure [71 we show the en¬ 
ergy conservation of the code for a given time-step accuracy. 
To avoid spurious energy errors upon start-up, we have mea¬ 
sured the energy errors in the interval from t = 4 to t = 5. 
The calculations become generally less accurate for larger 
time-step sizes as expected. However, especially for values 
of T-k ^ 0.3, the scaling is much steeper for larger values of 
?'chain,o- i-fi- when a larger fraction of tree code particles are 
integrated in the more accurate chain. For the simulations 
with Tchain.o = 0.02 and rchain.o = 0.03 the relative energy 
error scales roughly as oc t* . 


5 COMPARISON WITH ANALYTICAL 
ESTIMATES AND OTHER CODES 

5.1 Dynamical friction of a SMBH in a Hernquist 
sphere 

As a first test of rVINE, we investigate the orbital evolu¬ 
tion of a massive particle due to dynamical friction in a 
spherical non-rotating Hernquist sphere (model A), com¬ 
paring results from different simulations using NBODY7, 
rVINE, VINE, and Gadget-3. The massive particle of mass 
Mbh = 10 ~® Mtot is set on an initially circular, co-rotating 
orbit at the half-mass radius (rBH(t = 0) = 2.41) (see also 
Section (3]) . 


5.1.1 Dynamical friction theory 

For a meaningful assessment of the different codes, we com¬ 
pare the simulated SMBH trajectories wit h theoretical ex¬ 
pecta ti ons from dynamical f riction theory llChandrasekhail 
1 19431 : iBinnev fc Tremaiii^ l2008h . Assuming a locally 
isotropic velocity distribution function, a homogeneous den¬ 
sity p, as well as a sufficiently large velocity of the SMBH, 
ubh, relative to the stellar background, the rate of deceler¬ 
ation of the SMBH due to d ynamical friction may be g iven 
by the standard formula (cf. iBinnev fc Tremainell200^ Eq. 
8.5) 


ttDF = — 47r p(u < ubh) InA ubh, (6) 

■^BH 

where In A is the Goulomb logarithm. Only the mass density 
of stars moving slower than the SMBH, p(v < wbh), con¬ 
tributes to the dynamical friction. With gravitational soft¬ 
ening there is a limit on the maximum gravitational force 
that may be exerted in any star-SMBH encounter t hrough 
an effective minimum impact parameter bmin (see e.g. [White! 
ll976l:IJust et al.ll20ljl . 

bmin = 1.5 • max(£*, £bh). (7) 


Taking this softening effect into account, we write the 
Goulomb logarithm as 


In A = In 


\/bmin 


+ bnr 


( 8 ) 


where bmax is the maximum impact parameter and bgo is the 
impact parameter for a 90° scattering event of the incident 
star. The choice for the latter two parameters is often rather 
arbitrary, with bgo depending on the typical velocity utyp of 
the stars. The maximum impact parameter bmax is often 
taken pro portional to the orb i tal ra dius of the SMBH ran. 
Following ljust fc Penarrubial ll2005l l. we identify bmax with 
the local scale length, 


b 


max — 


p 

|Vp| 


r-BH 

3-p’ 


^ 2 , 


(9) 


where the last equation is true for t he family of p-models 
dPehnenl I 1993 I : llVemaine et al.l Il994[l if p ^ 2. Hence, in 
the case of the Hernquist profile (p = 2) we obtain that 
the maximum impact parameter exactly equals the orbital 
radius of the SMBH, i.e. bmax = fbh. In addition the 90° 
deflection parameter is given by 


bgo 


GMbh 




G Mbh 

2cr2 -b jj ■ 


( 10 ) 


Using Equations I® - ([I0|), we evolve the orbit of a 
SMBH, initially placed on a circular orbit at the half-mass 
radius, in the analytic Hernquist potential including the ad¬ 
ditional drag forces due to dynamical friction with a leapfrog 
integrator. The density of stars with velocities smaller than 
the SMBH is represented as 


p(u < UBh) = «s ■ ps{v < Ubh), (11) 


where pa(u < ubh) denotes a locally isotropic Maxwellian 
velocity distribution given by 


Ps(u < ubh) 


p(r) X 


rtv\ 

erf(A)-^e 


( 12 ) 


where X = UB H/v^a with dispersion g, and erf(A) is the 
error function l|Binnev fc Tremaine! 1200^ . We use Ks as a 
free parameter to account for the fact that the velocity dis¬ 
tribution of the Hernquist model does not follow a sim¬ 
ple Maxwellian distribution, as often used as an approx¬ 
imation in dynamical friction calculations. Hence, the lo¬ 
cally isotropic Maxwellian velocity distribution corresponds 
to Ks = 1.0. 


5.1.2 Results 

The upper panels of Figure [8] show the evolution of the 
radial decay (left) and the velocity (right) of the SMBH 
with time for model A_Nbody7_100k. NBODY7 is very 
well suited to follow the dynamical friction of the heavy 
body. We gauge run A_Nbody7 100k (blue line) against 
theoretical orbits for different values of Ks = (0.5, 0.75,1.0) 
shown as black lines and use it subsequently as a refer¬ 
ence for the other simulation codes. Initially the orbital 
evolution follows quite closely the analytic prediction with 
Ks = 0.75 (black dotted line). At smaller central distances 
(r < 0.6 t^b) the dynamical friction seem to act even more 
efficiently in the full N-body run. This leads to a slightly 
faster orbital evolution such that the SMBH reaches the 
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Figure 8. Comparison of the distance to the centre {left panels) and velocity {right panels) of a single SMBH, initially set on a circular 
orbit at the half-mass radius of a Hernquist sphere, for different code architectures. From top to bottom we show the direct summation 
code NBODY7 (blue; model A_Nbody7_100k), the AR-chain regularized tree code rVINE (red; model A_rVine_100k), and the VINE 
and Gadget-3 tree codes (green and pale blue lines) with varying softening lengths, £bh = 0.1 (models A_Vine-El_100k and A_Gadget- 
El_100k) and ^bh = 0.02 = e* (models A_Vine-E2_100k and A_Gadget-E2_100k). Theoretical expectations from dynamical friction 
theory are given as black lines for different parameters of k,s (/^s parametrizes the deviation from a locally isotropic Maxwellian velocity 
distribution, see Section fb.lll . For comparison, we show the results for the NBODY7 run also in the other panels. rVINE very well recovers 
the expected orbital evolution by effectively removing the limitations imposed by gravitational softening. Results shown are averages 
over several realisations of the initial conditions, as described in the text. 
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centre (defined as r < r{M < 2Mbh), where dynamical 
friction ceases to be efficient) within t < 500 time units. 
In the second row of Figure [8] we compare the efficiency of 
dynamical friction on the SMBH in model A_rVine_100k 
(red line) and model A_Nbody7_100k (blue line). We find 
that the early orbital evolution (r > O.Sr^jj) in model 
A_rVINE_100k is nearly indistinguishable from the one in 
model A_Nbody7_100k such that the SMBH reaches the 
centre on a very similar timescale, only about ~ 8% longer 
than in model A_Nbody7_100k, in very good agreement 
with the NBODY7 result. 

On the other hand, owing to their inherent resolution 
limitations, we expect the tree codes to show a significantly 
slower radial decay depending on the adopted gravitational 
softening, via Equation 0. If the SMBH is treated as a 
quasi collisionless particle with a large gravitational soften¬ 
ing length (ebh = 0.1, third row), it decays to only about 
half its original distance within the time span {t ~ 960) 
shown in the Figure |8] in both model A_Vine-El_100k 
(green line) and model A_Gadget-El_100k (pale blue line). 
The total time to reach the centre is t > 1350 time units 
for both codes, significantly longer (by a factor of > 2.7) 
than with NBODY7. This can be understood as a result of 
the reduced frictional force for near encounters. Both codes, 
however, follow the analytical prediction with Kg = 0.75 if 
we adopt the minimal impact parameter bmin = 1.5 • ebh 
from Equation 0. Setting £bh = e* = 0.02 (bottom row) 
yields an enhancement in the near encounters, but the decay 
time to the centre is still about 50% (70%) longer than in 
the NBODY7 case for model A_Vine-El_100k (A Gadget- 
E1 100k). Hence, even this drastic reduction in the gravita¬ 
tional softening of the SMBH does not significantly improve 
the accuracy of the dynamical friction time-scale. 

While all codes capture the effects of dynamical fric¬ 
tion as expected, if we consider softened gravitational force 
terms, this highlights the importance of an accurate treat¬ 
ment of close encounters for the correct description of dy¬ 
namical friction in N-body tree codes. From our analysis 
we conclude that the AR-regularized tree code removes the 
limitations due to the gravitational softening of the SMBH 
imposed on present-day tree codes by including nearby par¬ 
ticles in the chain. 


5.2 Evolution of a SMBH binary at the centre of 
a Hernquist sphere 

Another crucial test for rVINE is to investigate the harden¬ 
ing of a SMBH binary, compared to results obtained with 
NBODY7 and VINE. For our model B, we choose an ini¬ 
tial se t-up similar to the one presented in iMerritt et al.l 
ll2007tl and IVasiliev et al.l ll2014l l. In particular, we set up 
two SMBHs at the centre of a non-rotating Hernquist sphere 
with masses Mbh = 5 x 10“® Mtot and one of the SMBHs 
on a circular orbit (a;BHi = 0.1, Uy,BHi = Jucircl ~ 0.28) 
and the other SMBH at rest at the origirjf]. For models 

® Note that in IVasiliev et all (l2014h both SMBHs are set on co¬ 
rotating orbits, while here we set one of the SMBHs at rest at the 
origin. We tested that this different set-up gives similar results. 
Furthermore, in their simulations, the initial velocities are set to 
|u| = 0.31, which is about ~ 10% higher than the circular velocity. 


B_Vine and B_rVine we use gravitational softening lengths 
e* = 0.01 = ebh for the stellar and SMBH particles in the 
tree-code. Particles in the chain are not softened. The initial 
distance is chosen slightly larger than the influence radii of 
the SMBHs, which typically increase in the B_rVine and 
B_Nbody7 runs from ~ 0.1 to ~ 0.16 accord¬ 
ing to a decrease in the central density during the simu¬ 
lations. In the simulations using model B rVine, initially, 
only the co-rotating SMBH is regularized, but quickly - for 
t < 1.5 — the second SMBH is captured into the chain. In 
the B_Nbody7 runs we use an AR-chain with a chain ra¬ 
dius of Rchain = 1-25 X 10“®, while we test different values 
of the initial chain radius for the AR-chain in rVINE which 
we set to a size comparable to the hard binary semi-major 
axis, TchainO ~ & few X ahartd- The hard binary semi-major 
axis is given as 

= ^^ (13) 

Mbhi + Mbh2 4 

with fj, = (1/Mbhi + l/AfBH 2 )~^ being the reduced mass. 
Initially, Uhard = rf,^/16 Ri 6.25 x 10“^. 

In Figure [9] we show the time evolution of the SMBH 
binary hardness (l/fl, left panels) and eccentricity (right 
panels) for the first t = 100 time units in the simula¬ 
tions of model B with N = 10® particles. If we set a 
large initial chain radius, rchain.o = 0.02 (red solid lines), 
rchain.o = 0.03 (red dashed lines), or rchain.o = 0.04 (red 
dot-dashed lines), model B_rVine_100k agrees well with 
model B_Nbody7_100k (blue lines) showing an efficient 
hardening with a nearly constant hardening rate ^(1/a) for 
t > 10, and a rather mild binary eccentricity (e < 0.4). On 
the other hand, when setting the initial chain radius compa¬ 
rable to Uhard and the softening parameter, rchain.o = 0.01 
(red dotted lines), the hardening of the SMBH binary pro¬ 
ceeds much slower since we start to miss out on some close 
encounters with the hard binary, owing both to the effect of 
softening and the lower cross-section of the chain. Model 
B_Vine_100k (green lines), on the other hand, shows a 
qualitatively different picture, with a significantly larger sep¬ 
aration between the two SMBHs (1/a « 300) and a high 
binary eccentricity of e < 0.95 at the end of the simulations. 
While the reduced binary hardening is clearly due to both 
the softened gravitational forces near the SMBHs and the re¬ 
duced rate of loss-cone refilling in the collisionless tree-code, 
the reason for the high binary eccentricity is unclear. 

Figure [TO] extends the analysis of Figure to a range 
of different particle numbers. For VINE and rVINE we in¬ 
vestigate simulations with particle numbers ranging from 
V = 2 X 10'‘ - 10® and between V = 2 x 10'‘ - 10® for 
NBODY7. In the initial stages of the simulation {t < 5), 
when the loss-cone is full and dynamical friction is still effi¬ 
cient, the binary parameters evolve qualitatively similar in 
all three codes, with a steep rise in 1/a and a quite broad 
range of moderate eccentricities 0.1 ^5 e ^ 0.7). However, 
thereafter models B_Vine (bottom row) quickly evolve to 
high eccentricities while the binary semi-major axis typi¬ 
cally stalls at a < 0.5 tthard- In models B_Nbody7 (top row) 
and B_rVine (middle row), on the other hand, a hard bi¬ 
nary with a ^ Ohard forms quickly within t « 2.5 — 5. Models 

® Note that employing equally large chain radii would not be 
feasible in NBODY7 without major modifications to the code. 
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Figure 9. Evolution of binary parameters of a SMBH binary at the centre of a non-rotating Hernquist sphere (model B). Shown are 
binary hardness (1/a, left panel) and eccentricity {right panel) as a function time. Different colours indicate simulations using different 
codes, i.e. NBODY7 (blue), VINE (green, with sbh = £*)i and rVINE (red lines). Different line styles for model B_rVine denote different 
runs with varying initial chain radii (rchain.o = 0.01: dotted, 0.02: solid, 0.03: dashed, and, 0.04: dot-dashed). The results shown are 
averages over several realisations of the initial conditions, as described in the text. 


B_Nbody7 and B rVine again show an overall much more 
efficient binary hardening and eccentricity evolution than in 
B_Vine. The binary semi major axis and eccentricity reach 
final values of 1/a < (1100 — 2600) and e < (0.2 — 0.7) 
in the B_Nbody7 and B_rVine runs, respectively, depend¬ 
ing on N. Both in models B_Nbody7 and B_rVine the bi¬ 
nary hardening decreases with increasing N. Since this is 
due to the lower efficiency of collisional loss-cone refilling for 
larger N, however, the decrease is more pronounced in model 
B_Nbody7. As expected, in model B_Vine we again have 
a very weak evolution of the binary hardening {1/a < 400) 
but high final eccentricities (0.8 < e < 0.95) with a very 
weak dependence on N. 

In Figure [TT] we show the effect the ejection of low an¬ 
gular momentum stars at the bottom of the potential has on 
the structure of the galactic nucleus. Shown are radial pro¬ 
files for the density (upper panel), the radial velocity disper¬ 
sion (middle panel), and the velocity dispersion anisotropy 
parameter. 


Pa = l- 


I 2 

2 a/ 


(14) 


(bottom panel), for the models shown in Figure[9]at the end 
of the simulations. Stars are ejected from the galaxy centre 
by the hardening binary on orbits with high radial velocities 
in the B_rVine_100k (red) and B_Nbody7_100k (blue) 
runs. This leads to a large increase in the radial velocity 
dispersion at galacto-centric radii with r > lOro for these 
simulations (upper panel), together with a depression in the 
central density profile within r ^ and some added mass 
in the outskirts of the galaxy (middle panel). The central 
density profile is converted from an initial Hernquist pro¬ 
file with inner slope of ~ — 1 (dashed line) to a profile with 
p oc ® in the B rVine_100k and B_Nbody7_100k runs. 
The increase in the radial velocity dispersion is not seen in 


the tangential velocity dispersions such that the anisotropy 
profile is strongly radially biased for r > 10 ro. We verified 
that this is caused only by particles escaping the system 
after being ejected in interactions with the central binary 
in both codes. For B_Vine_100k (green), however, all ra¬ 
dial profiles evolves insignificantly throughout the simula¬ 
tion. rVINE seems on average to be more effective in re¬ 
moving mass from the centre than NBODY7. Given that the 
SMBH binary hardens by roughly the same amount in both 
codes, this might be due to the fact that the central density 
is replenished more effectively by the high rate of collisional 
loss-cone refilling with stars originating from larger radii in 
the B_Nbody7_100k runs. 


We further analyze the properties of the high radial ve¬ 
locity stars in the realisation using model B_rVine_100k 
shown in Figure [TT] by examining the radial velocities of 
stars escaping the chain in Figure [T2| Not all of the parti¬ 
cles interacting with the SMBH binary leave the chain on 
high radial velocity orbits: the majority of the escapers has 
Wrad < 1 at all times in the simulation. However, for t < 10 
there is an enhanced interaction rate of stars in the initially 
full loss-cone with the hardening SMBH binary leading to 
a significant population of escapers with radial velocities 
1.5 < Urad < 2.5, comparable to the expected kick veloci¬ 
ties in a slingshot interaction of a field star with a massive 
binary with hardness 1/a ~ 1000. Furthermore we find a 
small population of high-velocity outliers (2 < VraA < 6.5) - 
much higher than the expected kick velocity - which provide 
a clear observational signpost of the hard SMBH binary at 
the centre. 
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Figure 10. Time evolution of binary hardness (1/a, left panels) and eccentricity {right panels) as a function of particle number N for the 
three different codes in model B: NBODY7 {top row), rVINE {middle row) and VINE {bottom row). Shown are simulations with particle 
numbers increasing from N = 20k to N = IM from top to bottom, with colours indicated in the legend. The initial chain radius is set 
to rchain 0 = 0.02 in the simulations of model B_rVine. The results shown are averages over several realisations of the initial conditions, 
as described in the text. 
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Figure 11. Radial profiles of the final density {upper 
panel), radial velocity dispersion {middle panel) and anisotropy 
parameter {bottom panel) for the models shown in Fig¬ 
ure [9l B_Nbody7_100k (blue), B_Vine_100k (green) and 
B_rVine_100k (red lines). For comparison, the initial Hern- 
quist profiles are shown with the dashed lines, the dashed-dotted 
line shows an inner density profile following p oc and 

the vertical dashed lines indicates the resolution limit in model 
B_Vine_100k, given as ~ 2e*. Inserts show the residuals from 
the initial Hernquist profiles for the same radial range. Shown are 
averages over several realisations of the initial conditions. 


6 DISCUSSION 

Due to the favourable 0{N \og{N)) scaling of the tree code 
we should be able to employ (much) higher particle num¬ 
bers than in the test calculations presented in this paper in 
future applications of the new code. However, as a caveat, 
we note that the time spent for the AR-chain calculations 
scales steeper with particle number than the time spent for 
the tree calculations, mostly due to the costly extrapolations 
of the Bulirsch-Stoer method and the repeated predictions 
and force evaluations of the perturbers. For the simulations 
performed for this paper execution of the chain part of the 
code has only taken a moderate fraction of the total CPU 
time (typically below 5% with a maximum of ~ 20%), and 
with some optimisation it should be possible to further in¬ 
crease the speed of the code. There will, however, be some 
critical particle number, A^crit, at which the AR-chain will 
become the dominant contributor to the total computing 
costs even if only a small fraction of the total particle num¬ 
ber is actually integrated in the chain. It is not within the 
scope of this paper to investigate this in detail and we leave 
this to future work where we will use our new hybrid code 
for full-scale galaxy simulations. 

Throughout Section 15.21 we have found a qualitatively 
similar evolution of the hardening binary both in NBODY7 
and rVINE for chain sizes of ~ a few times the hard bi¬ 
nary distance flhard (Equation (I13II ). The agreement is par¬ 
ticularly good for the highest particle numbers studied here 
{N ^ 80fc) where spurious relaxation effects become less 
important in the direct N-body code. Similarly, for lower 
particle numbers {N 50k), rVINE shows a shallower N- 
dependence of the hardening rate since the tree-code better 
reproduces the collisionless galactic stellar dynamics at dis¬ 
tances far from the SMBHs. Hence, the hybrid code seems 
to catch the relevant dynamical interactions of a real galaxy 
better at lower N for our set of parameters adopted for the 
chain. 

However, we also note here that the high hardening 
rates we find in Section 15.21 in the NBODY7 simulations 
are somewhat larger than those found in a number of recent 
studies of spherical and axisymmetric m odels using similar 
techniques and initi al conditions (see e.g. iKhan et al.ll201^ ; 
IVasiliev et al.ll20l4l . In particular, the final binary hard¬ 
ness for our B_Nbody models are on average about a fac¬ 
tor of ~ 2.5 higher for compar able particle number s than 
the spherical models studied in IVasiliev et al. 1 ll2014h . Sev¬ 
eral reasons could be responsible for this difference including 
slightly different choices in the initial positions and veloci¬ 
ties of the SMBHs (see Section (5.211 or differences in the in¬ 
tegration techniques (parameters for the AR-chain, settings 
for the gravitational softening and the integration accuracy, 
etc.), the most likely being the different way of introducing 
the SMBHs in the initial conditions. This may lead to an 
overestimate of the hardening rate, especially at the begin¬ 
ning of our simulations. 

In Section ED we have found that in the tree codes 
VINE and Gadget-3 — even with the most conservative 
choice of the SMBH softening length — the dynamical fric¬ 
tion time-scales for the SMBH to sink to the galactic cen¬ 
tre differ by more than 50% from the ones in NBODY7 
and rVINE. Hence are present-day cosmological and galaxy 
merger simulations suffering from significant (unavoidable) 
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Figure 12. Radial velocity statistics of particles escaping the 
chain in one of the B_rVine_100k simulations, with ‘rchain,o = 
0.02. Top panel: radial velocities as a function of time, bottom 
panel: Distribution of radial velocities of particles leaving the 
chain, with the velocities binned by Av = 0.25. For c omparison we 
show the expected kick velocities (ukick oc p/a, ISaslaw et alj 
Il974h for inverse binary semi-major axes Ija = 100 (dashed black 
line), 1/a = 1000 (solid black line) and 1/a = 2000 (dotted 
black line), and the escape velocity from the centre of a Hern- 
quist sphere (dashed red line). 


uncertainties in the SMBH orbital time-scales? Strictly 
speaking they do, but probably, in most cases SMBHs are 
rarely found orbiting ’naked’ in their host galaxies, but are 
instead embedded in stellar and gaseous cores or cusps that 
are then the prime subjects to dynamical friction in galaxy 
interactions. However, it has been shown that ’naked’ BHs 
might be quite commonly formed after th e disruption of 
galac t ic nuclei in gas-rich minor mergers llCallegari et al.l 
I 2 OO 9 I : IVan Wassenhove et al.l l2014l ) , making accurate esti¬ 
mates of the dynamical friction time-scales of the ’naked’ 
BHs necessary in these cases. 

Of similar importance for the galaxy formation commu¬ 
nity is to get better estimates for SMBH binary coalescence 
time-scales in order to make robust predictions with respect 
to the dynamical evolution of SMBHs in their host galax¬ 
ies. For example, studying the exciting possible formation 
of systems with multiple SMBHs at high redshift due to 
the high merger rate of galaxies relies crucially on (1) an 
accurate description of dynamical friction in order to cor¬ 
rectly quantify the populations of binary, triple, etc. SMBHs 
being present at a given time, and (2) accurate orbits in 


order to reliably calculate the final outcome of the strong 
multi-body interactio ns between the SMBHs in the galac¬ 
tic c entres (see e.g. iKulkarni fc Loe9 120121 : iBlecha et al.l 
Obtaining accurate coalescence time-scales for binary 
SMBHs in gas-rich galaxy mergers is also essential for ac¬ 
curate estimates of the likelihood of recoiling SMBHs es¬ 
caping from the rapidl y steepening centra l potential of the 
merger remnants (e.g. ISiiacki et al.l 120111 ). This should be 
particularly relevant for large scale cosm ological simulations 
like e.g. the recent EAGL E and Illustris (ISchave et al.ll2015l : 
IVogelsberger et al.ll2014 1 simulations that assume fast coa¬ 
lescence of two SMBHs once their distanc e falls below the 
resolution limit (see e.g. ISiiacki et al.]l2014l . for details). 


7 CONCLUSIONS 

In this paper we have presented a hybrid code combining 
an OpenMP-parallel binary tree code (VINE) with an al¬ 
gorithmic chain regularization scheme, and report on first 
tests with the new code called “rVINE’Q. 

We have shown that, using the AR-regularized tree 
code, we can significantly improve the numerical accu¬ 
racy in the calculation of the gravitational interactions of 
SMBHs with their close environment. By comparison with 
the collision-less code NBODY7, we have verified that we 
have overcome some of the fundamental limitations imposed 
by the gravitational softening of the SMBHs, as it is used in 
traditional tree codes. As a consequence, we are now able to 
follow the orbital evolution of SMBHs much more accurately 
in more realistic, galaxy-scale settings. We have shown that 
with the new hybrid code we obtain both significantly im¬ 
proved estimates for dynamical friction time-scales of single 
SMBHs sinking to the galactic centres and for the time evo¬ 
lution of hard SMBH binaries. In particular, using rVINE, 
we find a clear N-dependence of the binary hardening rate, 
a low binary eccentricity along with a moderate eccentric¬ 
ity evolution, as well as the conversion of the galaxy’s inner 
density profile from a cusp to a core via the ejection of low- 
angular momentum stars on orbits with high radial velocity, 
similar to the results obtained with NBODY7 here and in 
previous work. 

Due to the modular design of rVine, the AR-chain part 
with its hybrid interface should be easily portable to other 
codes used for simulations of galaxy formation. It will like¬ 
wise be straightforward to incorporate additional physics 
into rVine, e.g. formulations for hydrodynamics and addi¬ 
tional sub-grid models of (1) star formation and stellar feed¬ 
back, and (2) black hole accretion and feedback, or, the ad¬ 
dition of post-Newtonian terms in the AR-chain. Note also 
that in the present paper we have restricted ourselves to the 
case of regularizing the dynamics of one single subsystem 
only. The next step here is to extend the present code to al¬ 
low for multiple chains in order to handle the regularization 
of several distant SMBHs at once. 

Important problems that will benefit from the accurate 
dynamical modelling of the evolution of (binary) SMBHs 
are predictions with regard to SMBH coalescence rates 


^ rVINE is available to anyone interested upon request to the 
authors. 
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and their asso ciated gravitational wave background (e.g. 
lHaehneitlll994ll . the population of SMBHs and AGNs (ei¬ 
ther living at their host galaxies centres or being displaced 
from the central regions by three-body encounters or grav- 
itatio nal wave re coils), the acceleration of hyper-velocity 
stars l|Hillsl fl988l l. and the imprint of SMBH binaries on 
the structural and kinematic proper ties of the inner parts 
of the stellar component of galaxies ( Ebisuzak^jtalJ Il99ll : 
iMilosavlievic &: Merrittll200lK iMeiron fc Laoil 2013t l. 

Further improved hybrid codes like the one presented 
here will help to bridge the gap between different fields of 
astrophysics that are currently still separated by huge dif¬ 
ferences in the relevant scales of space, time, and mass and 
should eventually allow to study properly the effects of stel¬ 
lar dynamics in the sphere of influence of central supermas- 
sive black holes on the structure of their host galaxies. 
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